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' The recently introduced models of reionization bubbles based on 

' extended Press-Schechter theory (Furlanetto, Zaldarriaga & Hern- 

, quist (2004a)) are generalized to include mergers of ionization sources. 

\^ I Sources with a recent major merger are taken to have enhanced photon 

■ production due to star formation, and accretion onto a central black 

c~| . hole if a black hole is present. This produces a scatter in the number of 

Oh| ionized photons corresponding to a halo of a given mass and a change 

O ' in photon production over time for any given halo mass. By extending 

+l3 . previous methods, photon production histories, bubble distributions, 

^ I and ionization histories are computed for several different parameter 
and recombination assumptions. The resulting distributions interpo- 
late between previously calculated limiting cases. 



1 Introduction and Background 

Reionization marks the historical event when hydrogen in the universe trans- 
formed from mostly neutral to mostly ionized; it has only recently become 
accessible via observations and numerical simulations and has stirred great 
interest and advances in our unde rstand ing of early struct ure for mation (see 



e.g. revi ews b y Barka na fc Loeb (|200l[ ). Loeb & Barkana (|200l[ ). Cooray & 
Barton goflS), Loeb goQ^). 



Current observations suggest that the process of reionization is complex, 
perhaps lasting over an extended period of time (for a review, see Fan, Carilli 
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& Keating tOO(t). The optical depth to electron scattering seen by WMAP 



(Spergel et al ()2o'o6[ l ) can be modeled with reionization starting by z > 10. 
On the other hand, spectroscopic observations of high-redshift quasars find 
Gunn-Peterson absorption due to int ervening neutral hydrogen in the inter- 



galactic medium (Becker et al (j2nnih . Fan et al (2003), White et al ^005)), 



although w ith si gnificant fiuctuations along different lines of sight (Oh & 
Furlanetto ( 2005[ )). which suggests that the tail end of reionization occurs at 



z ^ 6. 

Theoretically, questions about reionization range from the fundamental 
timeline of reionization to detailed characteristics. The latter include the 
nature of ionizing sources, the topology of the ionized regions and their 
evolution, and the behavior of the ionized regions relati ve to matter over- 



densit ies an d voids. Numerical approa ches ( e.g. Gnedin ()2000), Ra zoumov 
et al (j2002h . Ciardi, Stoehr & White (j2003h . Sokasian et al (|2003l : Eo04l ) 



are fruitful due to the complex nature of the sources and nonlinear cluster- 
ing involved. For small regions, simulations can track the nonlinear effects, 
although much known physics is too difficult to incorporate, and much un- 
known physics remains. However, the large disparity in size between the 
ionized regions (which can grow to be tens of Mpc comoving, e.g. Wyithe 
& Loeb (2QM), Furlanetto, Zaldarriaga & Hernquist ( 2004a| )) and the small 



scale distributions and physics of sources producing the photons is difficult 
to capture in a simulation. Recent numerical simulations have just begun 
to combine large and small scale effects in lar ger vo lumes (see for instance 
Kohler, Gnedin & Hamilton (j2005|), Iliev et al (1200,51 )). 



Analytic models are useful to encompass the wide range of scales in or- 
der to try to identify generic behavior. In addition, analytic models allow 
exploring a larger range of parameters and other assumptions. This paper 
extends an analytic model introduced by Furlanetto, Zaldariagga and Hern- 
quist (|2004^ (hereafter FZH) to characterize the growth of ionized regions 



or bubbles. The basic idea is to move beyond the initial step of cons iderin g 
HII regions aroun d indi vidual collapsed halos (e.g. Arons & Wingert 



Barkana & Loeb (|200lh to consider regions around several collapsed halos. 



The simplest form goes as follows: collapsed halos of mass m are able to 
ionize a mass M where 

M = Cm (1) 

and C is a constant. Halos need to be massive enough to produce photons, 
thus m > mmin{z) = m(T = 10"^ K, z) is imposed, i.e. they need to be massive 
enough to have efficient atomic Hydrogen line cooling. For a constant C,, 
independent of m, the size of ionized regions can be calculated immediately 
in closed form: an ionized region of mass M has the fraction of mass fcoii 
bound in halos above mjnin{z) equal to or greater than Using extended 
Press-Schechter theory (references below), the mass fraction fcou in a fully 
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ionized bubble of mass M then obeys 

hou = erfc[ ^ ^ - ] > . (2) 

This constrains the average overdensity 5m for any M . Here Cmin denotes the 
density fluctuations smoothed on a scale of IXlclSS '^^rnini^^ 

), likewise for (j^(M), 

and 5c{t) is the threshhold density for collapse. The overdensity required for 
collapse is 5c = 1.686, in this picture the densities stay constant and the 
threshold lowers with time, so 5c{t) = 5c/D{z{t)) where D{z) is the growth 
factor. The bubble's average overdensity relative to 5c{t) is 5m, corresponding 
to physical overdensity 5R^M{t) = 5MD{z(t)), and Vm(1 + 5R^M{t))p = M, 
where Vm is the bubble volume and p is the mean density. The condition for 
a bubble to be ionized, equation |21 can be rewritten in a form that is more 
easily generalizable, 

Cfcoll = 1 . (3) 

The function 5m allows one to determine the number distribution of bubbles 
of mass M and other properties discussed below. 

In FZH, recombinations were included implicitly in (. A more sophisti- 
cated approach balancing the total number of ionizations and recombinations 
per unit tim e in t he intergalactic medium (IGM) was introduced in Furlan- 
etto & Oh ( 2005[ ) (hereafter FO05), where they considered two models of 



IGM gas density distribution for recombinations: one model adopted an an- 
alytical fit of gas densit y to simulations at low redshift by Miralda-Escude, 
Haehnelt & Rees and extrapolated it to high z; a second model con- 



sidered minihalos as photon sinks. Minihalos are neutral, dense blobs with 
mass IO^Mq, randomly distributed in the IGM. 

These results were extended by Furlanetto, McQuinn & Hernquist ( 2005[ l 



(hereafter FMH05) and others. For example, FMH05 include the effects of 
different halo mass functions, a mass dependent C("^)) cind stochastic fluc- 
tuations in the galaxy distribution. They found the consequences of these 
additional effects on the bubble size distributions and corresponding evolu- 
tion in time, the emissivity inside a bubble, the observable neutral hydrogen 
21 cm power spectra, and the kinetic Sunyaev-Zel'dov ich sig nal. Tools for 
calcul ations of the latter were developed by Zahn et al (20051) and McQuinn 



et al (|2005a[) . they also extended many aspects of the approach. Other con- 
sequences of this model for observables have b een calculated, for example 
in Furlane tto, Hernquis t & Zaldarriaga (12004 ) . Furlanetto, Zalda rriaga & 



Hernquist (|2004bl : l2006h . McQuinn et al (|2005bl ) and Alvarez et al (|200.5l ). 



Here we extend this model further to include the effects of major merg- 
ers. The original models are based on a time independent ( which only 
depends upon mass. Halos which have recently had a major merger produce 
an increased number of photons, due to starbursts (1:3 halo mass ratios) or 
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induced black hole accretion (1:10 halo mass ratios). Consequently major 
mergers introduce a time dependence in the average of ( and a scatter in (. 
We expect halo major mergers to be important, as even the smallest objects 
producing photons (chosen to be those with T > lO'^K) have masses several 
orders of magnitude above M^, at high redshifts, and thus are highly biased 
and very actively merging.^ We use estimates of photon production due to 
quiescent star formation, merger induced starbursts and black hole accretion. 
As the photon production estimates have many unknowns, we are especially 
interested in features which characterize the merging itself rather than the 
details we choose.^ 

There are several related studies in this active area of research. This 
is not the mergers of the bubbles, which was already co nside red in FO05. 
The scatter discussed is not the scatter (Furlanetto et al ()2005| )) in C due t o 



the stochasticity of the number of collapsed halos (Barkana & Loeb (j2004a|), 
Babich & Loeb (,2005|)), or the d ifference between metal rich and metal poor 
regions (Barkana & Loeb ^200^ ). Other related work includes the calcula- 



tion of the effects of mergers in producing cumulative numbe r of photons by 
qu asars up to redshifts ~ 5, e.g. by Wyithe & Loeb and Madau et 



al (|2004( ). Here we want the effects of the mergers in producing a time de- 
pendent photon production rate, with scatter, for collapsed halos of a given 
mass, and the consequences for the enclosing bubbles. 

Section §2 gives the analytic framework and formulae for the merger rates, 
recombinations and scatter. The natural quantity to use is the instantaneous 
photon production rate which then determines the time derivative of (t rather 
than ( itself. The resulting formalism is a natural extension of earlier work 
(e.g. FZH04, FO05, FMH05): their early time limit took reionization to only 
depend on the halos present at the time under consideration and included 
recombinations implicitly, the late time limit imposed equilibrium between 
instantaneous photon production and recombination rates. The extension 
of these earlier formalisms, introduced here, finds the total number of ions 
present by integrating the photon production rate over time (due to the halos 
present at any instant), and subtracting recombinations. Three prescriptions 
for recombinations are used. The scatter requires additional assumptions, 
we describe and choose a simple case next. In section §3, we calculate the 
analytic form of the time derivative (t in terms of the star formation and 
black hole accretion parameters and then explicitly show its dependence upon 



^M^t, a function of time, is the mass at which a{M) — 6c{t). An estimate of the 
formation minus destruction rate for ha los of mass M is proportional to 6c{t)/a'^{M) — 
l/dc{t), see e.g. Kitayama & Suto l)l996j) . For M > this rapidly increases with M and 
can be thought of as indicative of not only a high formation rate but also a high major 
merger rate-more detailed calculations confirm this trend. 

^Ideally this could eventually be turned around to constrain assumptions going into 
the photon production estimates but it is not clear how constraining that might be given 
the vast uncertainties. 
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mass and redshift for a few examples. Section §4 shows the effects on bubbles 
using a fiducial model and some variations: the bubble overdensity 5m as a 
function of bubble mass M (as a function of redshift and initial conditions), 
the characteristic bubble sizes and ionization fractions over time Xj, and the 
bubble size distributions associated with the (major) merger scatter. Section 
§5 concludes. An appendix spells out the basic Extended Press- Schechter 
quantities, details the estimates used to find the minimum halo mass for a 
black hole to be present and describes some other possible assumptions for 
scatter. The term ionized mass fraction refers to the mass found in ionized 
bubbles unless otherwise stated. 



2 Analytic Framework 

In this section we describe the analytic calculations, assumptions and formu- 
lae used in the rest of the paper: the halo mass function, major merger rates, 
recombinations and the merger related scatter. The input parameters to most 
of our calculations are VLm = 0.3, VL^ = 0.7, erg = 0.9, flbh'^ = 0.02, h = 0.7 a 



scale invariant initial power spectrum n = 1, and the Eisenstein-Hu ()l997l ) 
transfer function. Our fiducial model has these parameters. We considered 
several changes of cosmology. We used the above and only changed n to 1.05, 
which enhances fluctuations for large k, i.e. small scales, and we fixed to the 
above but instead took Qbh"^ = 0.225. Another variant we considered was 
the set of parameters from the recent WMAP analysis (Spergel et al (j^06)): 
= 0.24, = 0.7, as = 0.74, n^h'^ = 0.0223, h = 0.73, we discuss these at 
the end. Masses are defined in terms of h~^MQ and all lengths are comoving. 

2.1 Halo Mass Function 

The starting point for the extended Press-Schechter bubble model of reion- 
ization is the number density of col lapsed halos of dark matter. We use 
the Press-Schechter formalism to estimate numbers and source ma- 

jor merger rates. Examples of the resulting number densities and properties 
have been detailed by Barkana & Loeb (2001) and Mo & White (12002). 

Our use of the Press-S chechter mass function warrants some discussion. 
The Sheth-Tormen (|2002h mass function is a better fit to simulations at low 



redshifts; mas s fun ctions found i n high redshift simulati ons (J ang-Condell 
fc Her nquist (j200ll ). Reed et al (|2003t ). Heitmann et al (|2006l ). Ihev et al 



( 2005| )) often lie between Sheth-Tormen and Press-Schechter (the bes t agre e 



ment for Heitmann et al was with the fitting function of Warren et al ( 2005[ l). 



Agreement with dark matter simulations is not necessarily indicative of ap- 
propriateness either, as baryons and da rk ma tter are not as tig htly c oupled 
at high redshifts (e.g. Naoz & Barkana ^QQ^, Bagla & Prayad (|2006l) ). For 



6 



the original model (which we modify here), the the effect of using the Sheth- 
Tormen mass function instead of the Press-Schechter m ass fu nction for the 
number of sources was seen to be small (Furlanetto et al ( 2nn5h ). We altered 
the "tilt" n in one set of runs from n = 1 to n = 1.05 to see the effect of 
increased small scale structure as found in these simulations. 

The Press-Schechter mass function is the most tractable and our interest 
is in trends rather than exact numbers. With it, the number of recently 
merged halos can be calculated easi ly using exte nded Pr ess Sc hechter theory 
(Bond et al (|l99l[) . Lacey & Cole (|l993l : Il994[ ). Bower (|l99lh . Kitayama & 
Suto ()l996t )).^ There are problems as well with extended Press-Schechter 
theory, especi ally w hen considering large mass ratios in merger rates (see 
Benson et al ()2005| ) for a recent discussion of this). We restrict ourselves to 
major mergers, with small mass ratios. In addition to being in the regime 
where extended Press-Schechter works the best, this allows predecessors and 
final halos to be identified uniquely (a halo only has one predecessor with 
mass greater than half the final halo mass). 



2.2 Photon Sources 
2.2.1 Major Merger Rates 

After a major merger, a halo produces extra photons for some specified 
amount of time, denoted tion- We are interested in the number of halos 
of a given mass which are still "active," i.e. which are still producing extra 
photons due to a recent major merger. We ignore the lag between the ac- 
tual merger and the beginning of the photon production, which should not 
change the general trends of interest. We also assume that the extra photon 
production rate is constant during tion after the merger, so that we do not 
care about when exactly the merger happened, only that it happened within 
this relaxation time. 

The number of recently merged halos with mass m at time t is the product 
of two factors, depending upon two times: the time of the merger tj and the 
subsequent time of observation t, where the latter is close enough to U for 
the halo to still be excited from the merger, i.e. t — ti< tion- The first factor 
is the fraction of halos that have mass mg which have jumped at time tj from 
mass mj, Pi(mj — >■ mQ]ti)^dmidti. (The overdot denotes derivative with 
respect to t,.) The factor of mo/rrii makes the conversion from the number 
of points coming from rrii halos (the quantity given by the formalism) to 
the total number of points in their descendant mo halos. ^ We then want 

■^These are actually formulae for fast mass gain and can include accretion. Analytic 
for mulae differentiating betwe en ma jor mergers and ac cretion have been found by Raig et 
al l|l998l) . Salvador-Sole et al l|l998l) and Andreu et al l|2n0lh . 

^If the mo halo has come from an rm halo, all of its mass was not previously in the rrii 
halo, only rrii /mo of its mass. But the quantity wanted here is the full mass of halos which 
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this "excited" halo of mass mo to be in a halo of mass m at time t. The 
probability that a point in a halo of mass m at time t was in a halo of mass 
mo at time ti is Pi (mo, ti|m, t)dmo- Formulae for Pi and its derivative are in 
the appendix. 

For bubbles we need to go one step further and consider the regions of 
mass M surrounding these mass m sources. The probability that a halo 
(bubble) of mass M and ovcrdcnsity 5m contains a halo of mass m at time t 
is denoted by Pi(m, t\M, 6M)dm, corresponding to a number density of halos 
n{m,t\M,dM)dm = ■^Pi{m,t\M,dM)dm. (When a time t rather than 6 is 
written as an argument of Pi the implicit assumption is that one uses 6c{t) in 
Pi, the threshold density corresponding to time t.) Thus the number of halos 
which have merged at time ti from mass m^ to mass mo and are in a halo of 
mass m at time t within a bigger halo (bubble) of mass M and overdensity 
Sm is given by (see Appendix for more discussion and details): 

VMn(m, t\M, 5M)Pi{nT'i — ?7io, ti)Pi(mo, ti\m, t) — dmdtidmidmo ■ (4) 

mj 

2.2.2 Total Photon Rates 

The rate photons are produced in a region with mass M and overdensity Sm 
is the number of additional photons due to the mergers, plus the quiescent 
number of photons due to the collapsed halos present even if no merger 
has occurred. (We assume that each photon produced corresponds to an 
ionization in the region, in principle the bubble region can have some "escape 
fraction" which would effectively lower the number of ionizations in a region 
relative to the number of photons produced.) A constant C for sources of 
lifetime AT gives a Ct = C/^T. More generally, there is a dependence 
upon the mass of the initial and final masses before the merger [mi, mo), 
the mass m at the time of interest, and ti,t to give (t{nT-i,n"i-o,ti,m,t). The 
quiescent rate is taken to be Ct,q{^)- These quiescent photons are from stars 
not created during starbursts, i.e. ongoing star formation, in part due to 
accretion. Putting this together to get the rate mass is ionized gives 

/ dmn{m,t\M, 5m) j 

{/ dtidmidmoCtirUi, mo,ti, m, t)Pi{mi mo, ti)Pi{mo, ti\m,t)^^+ Ct,g(H} 
^ J dmn{m,t\M,5M)jCt,a{'^) 
= (C/)t 

(5) 

The above defines {Cf)t and C,t,a{'m)- Section §3 derives estimates for 
C,t{mi,mo,ti,m,t) using models for star formation and black hole accretion, 

have had mergers, i.e. halos which now have mass mo but had such an m, component 
earher; thus the factor mo/mj is included. 
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and gives parameters and limits for the integrals. In the integrals, the mass 
limits are determined by the mass ratio criteria for major mergers. In all 
cases we require mi^max/fno > rrii/mo > 0.5. The lower limit means some 
major mergers are neglected, but ensures that no major merger is counted 
twice, as only one halo with rrii > mo/2 can end up in mo. This is thus a 
lower limit on the number of mergers. The time limits in the integrals are 
determined by the time scales of relaxation after the mergers. Both the mass 
and time limits can differ for starbursts and black hole accretion. 

In principle mo,mm can be as small as m(T = 10'^K,z) = "^^^^^-^1% H'^Mq. 
However the probabilities above double count a halo which has two mergers 
within a relaxation time, which is more likely if mo,mm is chosen to be small. 
For the number of photons contributed, this is perhaps accurate (but perhaps 
not, e.g. for starbursts there might be gas missing after the first merger for 
some period of time). However for the scatter it will give two final halos with 
recent mergers rather than just the one which had two mergers. However, if 
'^o.mm is taken to be too large, the number of photons will be undercounted 
and the scatter underestimated. We used both mo,mm = 0.1m and mo,m.m = 
mminiT = lO'^ K,z) in our calculations,^ for nontrivial mass dependence of 
star formation (referred to as a = 2/3 below) these two choices are essentially 
indistinguishable. 

The ionized mass fraction in a region of mass M with overdensity 6m, the 
generalization of equation 01 is then 



where Rrecomb{t) is the instantaneous recombination rate in the bubble vol- 
ume, depending upon Rionit), the effective radius of the ionized region at 
this time. The dependence upon Rion takes into account that not all regions 
inside a bubble at time t are ionized at earlier times: recombinations can 
only occur where there are ions already present. As the first term gives the 
total number of photons produced in the region, assumed to produce the 
same number of ions, recombinations can be directly subtracted off. 

^If we take the limits of toq to be independent of m, then the integral over m can 
be done explicitly using an identity for composing the probabilities. This gives the mass 
fraction which had mergers at ti to mass tuq in a halo of mass M as 





/ 



dmPi(mi — !■ TriQ'jti) — Pi{rno,ti\m,t)Pi{m,t\M, 6]\i)d'mid'mQdti 



(6) 



III 

= Pi{mi —> mo;ti) — Pi(mo,ii|M, 5M)dm^dmodti 



(7) 



(8) 



i.e. it doesn't matter what m halo the original toq halo is in at time t, just that it is in 
the bubble M. 
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2.3 Recombination 

Recombinations decrease the number of ions present. We use the two esti- 
mates for the recombination rate for an ionized region of radius Rinn ( t) and 



overdensity Sm chosen by FO05 and a third by Melhma et al ( 2006| ). The 
recombination rate per hydrogen atom inside an ionized region of radius Rion 
at time t and overdensity 6m is 

aAiT)neCiRionit)) = + 5ij,A/(t))C(i?,o„(t)) , (10) 

where ne = ne(l + 5/?,m(^))- The physical overdensity Sji^Mit) enhances the 
rate at mean density = Q;^(T)ne = 2.4 x 10^/Myr. (We use a^(T = 
10'^ i^) for case A recombination, case B recombination (optically thick) takes 
Au — >■ 0.6^4^.) The three recombination models have different clumping 
factors, described in subsections below. 

To go from the ionization rates per hydrogen atom in the bubble to the 
change in ionization fraction in equationElrequires the comparison of counts^: 
counts of ions present and counts of recombinations depleting ions. The 
ionized mass fraction times the hydrogen mass of the bubble gives the number 
of ions in the bubble at a given time. The recombination rate at the same 
time is the rate per hydrogen atom times the number of ionized hydrogen 
atoms, nh{l + 5i?,A/(^))^on- Here Vion is the volume of the region where the 
hydrogen is ionized and Mion is its mass. Dividing by the total mass of the 
bubble and noting Vion = "^^aP {i+5R M{t))p ^ resulting ionization fraction 
with recombinations included, at some time t, is 

dt'[iCf)t' - Au^^il + SRMt'))CiRionim . (11) 

The mechanics of doing this integral are discussed in section §4, we now 
describe the clumping factors C{Rion) for the three different recombination 
models. 

2.3.1 The MHR model 

The first model (hereafter called the MHR mod el) is a smooth IGM gas 
distribution of Miralda-Escude, Haehnelt & Rees ()2000| ). They use the vol- 
ume density distribution of IGM gas, -Pv'(A) (A = p/p) fit to simulations at 
z ~ 2 — 4, finding 

/A -2/3 _ f )2 

Py{A)dA = AoA-^ exp[- ^ 2(2(5o/3)^ ^^^^ 

with jS = 2.5. Aq and Cq are set by requiring mass and volume normalization 
and 6q is the variance of density fluctuations smoothed on the Jeans scale for 

^We thank S. Furlanetto for discussions about this. 
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an ionized medium (at higher 60 = 7.61/(1 + z) to better than 1%). The 
distribution Py(A) is taken to be independent of environment (rather than 
depending upon 6m)- For recombination, they assume all gas below some 
density threshold A < Aj is ionized and everything above is shielded. One 
finds Aj by noting that recombination limits bubble growth- i.e. the mean 
free path Aj = Ao[l — Fy(Aj)]^^/'^ is the radius of the ionized region Rion- 
Here XoH{z) = 60km s^^ (in physical units) and Fy(Aj) is the fraction of 
volume with A < Aj. (If the region's radius Rion < Aq then the assumption 
is that recombination is negligible and Aj is set to zero.) The factor C is 
then 

CMHR{R^on) = dAPv{A)A^ , (13) 

JO 

again, note R^on enters implicitly via Ai(i?io„). At large Rion CMHR{Rion) 
tends to asymptote to a constant, and for most (all) radii it is smaller than 
Cmh{Rion) {Cm IPs) described below. 



2.3.2 Minihalos 

The second model is a minihalo model (FO05) of small dense absorbing 
clumps, taken to have mass of Mn,h = 10^ Mq with comoving mean free path 



- 15.7M^)V3(MJ)(^)2/3(^y!)-i/3h~i Mpc 



(14) 



and 

Cmh{Rion) = (1 — fmh) GW{Rion/ ^mh) (15) 

where fmh = 0.05 is the mass fraction taken to be in minihalos. Unlike the 
other two models it is redshift independent. Note that the effect of minihalos 
on recombination and thus reionization is likely to be more complicated than 
this simple prescription, which is pr obably an overestimate (Iliev, Scanna- 
pieco & Shapiro (,2005^ . Ciardi et al (|2006t) ). 



2.3.3 The MIPS model 



A third prescription from clum ping, m atched to numerical simulations is by 
Mellema, Iliev, Pen & Shapiro ()2006h (MIPS), which is an improved fit for 



the IGM clumping factor to the one given by Iliev et al ()2005l ). They used 
a A^-body simulation with a 3.5 Mpc box and 3248^ particles to get a 
mass resolution down to the Jeans mass; they then took out all found halos, 
including minihalos, to eliminate their contributions to the IGM density field. 
The resulting IGM was fit to give an Rion independent clumping factor: 

Cmips = 27.466 exp(-0. 1142 + 0.001328^^) (16) 
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It resembles the MHR clumping factor at large Rion as the former asymptotes 
to a constant for large Rion and high z, but is larger numerically, and thus 
more effective at suppressing bubble formation than the MHR model. This 
will be referred to as the MIPS model in the following. 



2.4 Scatter 

Extended Press-Schechter gives the fraction of halos which have had recent 
mergers. This is an average quantity, just as the number of collapsed halos 
in a larger region of mass M is given only on the average by n{m, t\M, 6m)- 
In any region a scatter in the number of sources will lead to a scatter the 
number of ionizing photons and thus in bubble sizes. For n{m,t\M,6M), 
FMH05 took the scatter due to the stochastic nature of the collapsed halo 
distribution (halos of mass m) in the larger bubbles of mass M to be Poisson. 
Numerical simulations have found this scatter to b e with in a factor of two of 
Poiss on for the cases studied by Sheth & Lemson (199^ and Casa-Miranda 



et al (j2nn2h . In equations, the ionized mass fraction is 

/in 
dm((m)—n(m,t\M,SM) (17) 
P 

and they assumed for the scatter 

{n{m,t\M,5M)n{m',t\M,5M)) = n{m,t\M, 5M)n{m' ,t\M, 5m) 

+ '-''^^nim,t\M,SM). ^ ^ 

where 6d is a Dirac delta function. 
For our case we want the scatter of 

/TTl 
dmn(m,t\M,6M)—Ct aim) ; (19) 
P ' 



I.e. 



TTlTn 

dtdt'dmdm'——{n{m,t\M,6M)Cta{m)n{m',t'\M,6M)Ct'a{m')) ; (20) 
P 

Thus instead of the scatter of n{m, t\M, 6m) for fixed time above, we need 

(n(m, t\M, 6M)n{m', t\M, 6m)) ■ (21) 

In addition, Ct,a{m) has scatter. To our knowledge the appropriate distribu- 
tions have not been calculated in numerical simulations, and thus additional 
assumptions are required to continue. 

The simplest assumption is to take the distribution of n(m, t|M, 5^/) to 
be Poisson both in mass and in time: 

(n(m, t\M, 6M)n{m', t'\M, 6m)) = n(m, t\M, 6M)n{m', t'\M, 6m) 

+^6nim - m')6D{6{t) - 6{t'))n{m, t\M, 6m] 

(22) 
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Note that n is a function of t only through S{t) and thus S{t) is taken to be 
the argument of the Dirac delta function. Using doit — t') would give the 
wrong dimensions. 

The second assumption is for the scatter of 

/TTl 
dtidmidmoCt{rni,mo,ti,m,t)Pi{'mi mo,ti)Pi{mo,ti\m,t) Kt,q{m) . 

723) 

We define ?T,°(mo|m, t) which appears in Ct,a(^) the following way: 

dtidmiCt{mi,mo,ti,m,t)Pi{mi ^ mo^U) — Pi(mo, ti|m, t) = V — n''(mo|m, t)C"(mo) . 

a=l P 

(24) 

The number density n°'{mQ\m,t) counts the recently merged halos of mass 
mo in a halo of mass m, with the index a denoting the range of integration 
for mi,ti. There are three ranges, corresponding to starbursts, black hole 
accretion or both. For each a, Ct(mj, mg, tj, m, t) only depends on which of 
these three cases is at hand and the values of mo,m, and so it is denoted 
C"(mo). Values for mQimo) are derived in the section and displayed in 
Table 1 for the starburst and black hole accretion case (the sum of these is 
needed when both are present). The Poisson assumption for n"(mo|m, t) is 
then 

(n"(mo|m, t)n''(m[,|m, t)) = ~ """^ ^"V(mo|m, t) (25) 

There is also an assumption being made that the number with scatter is 
related to the product of Pi Pi, rather than a separate scatter for the merger 
rate to mass mo and then scatter for the inclusion of these mo halos in m. 
The full scatter is then 

(G,a(m)n(m, t\M, SM)Ct'A^'Mm', t'\M, 6m)) 
= Ct,a{^)^i^^ ^1^' 5M)Ct',a{m')n{m\ t'|M, 5m) 
+^5B(m - m')5nm ' <5(t'))n(m, t|M, 5m)(CLM) 

= 6,a("^)^("^, t\M, SM)Ct'A^')^i^'^ ^'1^5 ^m) 

+^SD{m - m')SD{S{t) - 6{f))n{m,t\M,SM){Clai^) 

+ J dmodmidti{l + 5R{t,ti)){C,t{mi,mQ,ti,m,t)m)'^-^Pi{mi mo,ti)Pi{mo,ti\m,t)} 

(26) 

Here m = (1 + Suit, ti))pVm was used, where 5_R(t, ti) is the physical overden- 
sity 6{t)D{z(ti)) and the combination (t{^i,mo,ti,m,t)m was pulled out. 
We then define the scatter in / dt{(f)t via 

A( J iC f),dt)' = { J dtdt\Cf)t{Cf)t') - { J dt{Cf),){ J dt\a)t') (27) 

We calculated the consequences of this simplest set of assumptions for our 
models; some other possibilities are described in the appendix. 
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The scatter is an integral over time. In practice we cannot integrate back 
to arbitrarily early times, and so an initial condition and initial scatter are 
needed. See below for discussion. 



3 Photon increase from mergers 

The analytic expressions above require estimates of Ct{mi, mo, ti, m, t) and the 
quiescent rate Ct,q{Tn) to give concrete results. The two main sources of extra 
photons due to a major merger are starbursts and, if a central black hole i s 
present, accretion onto a central black hole (Kauffmann & Haehnelt 



Cavaliere & Vittorini (j2OOO0). Merger induced star formation and black hole 



accretion at high redshift might differ from their low redshift counterparts, 
detailed calculations and simulations have not yet been done. We use low 
redshift calculations and measurements as a guide. We stress that we are 
most interested in the effects of the time dependence, as the unknowns of the 
detailed modeling are vast. 

A major merger to a halo of mass m will add a total number of photons 
N^{m) over a time tjon- These photons can ionize the hydrogen in a 
region with mass M obeying i\L = J^^lh2i where X = 0.76 and m„ is the 

' Tip i Ifji ^ 

proton mass.'' This is the amount of mass M that can be ionized over the 
whole time T, giving, if M = (m, 

C = (28) 

We will want to compare the contributions to photons at some given time, 
thus we will look at all the photons being contributed assuming a constant^ 
production rate, i.e. (t = -r~- The times tjon for black hole accretion and 
starbursts differ in general. 

3.1 Star Formation and Starbursts 

Photon production du e to s tar formation can be estimated via (e.g. Loeb, 
Barkana & Hernquist ( 2nn5[ )) 



dN^, stars _ m Qb A^^on ^29) 



^In practice there will be a small correction, perhaps of order 10% due to the possibility 
that some of the photons will ionize Helium instead, which should be kept in mind. It can 
be thought of as a rescaling of our fiducial input parameters. We thank the referee for 
pointing this out. 

^More accurately, the accretion onto a black hole, if Eddington as is commonly assumed, 
tends to increase with time until all the infalling material is gone, while the starbursts 
tend to decay with time until all the gas is gone. However other estimates of black hole 
accretion have the accretion also slow down with time as the gas supply is decreased. 
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leading to 

Ct{mi, mo, ti, m, t) = — — ^ (stars) (30) 

Here Nion is the overall number of ionizing photons per baryon in the halo. 
We assume that the starbursts change the rate from a constant quiescent 
rate to a constant starburst rate for a time + thurst, until the stars due to 
the starburst are gone.^ 

For Nion our parameters are similar to Loeb et al ( 2005[ l. We take a star 



lifetime t^, = 80Myr. Showing our choices in parentheses^" , the efficiency 
with which baryons are incorporated into stars is /*(=0.1), and each baryon 
in a star produces N^^bary (=4=4:00) photons, with /esc(=0.05) of them escaping. 
Combining these we get 

Nion N^ i,(iryf*fesc ^ ^ ^ref(z) (31) 

The power a = 2/3 is to take into account the scaling of star formation found 
in low stellar mass galaxies (below msteiiar = 3/i x lO^^h'^Mo) by Kauffmann 
et al 6003), we assume that it is the same for both starbursts and quiescent 
star formation as in the z = sample low mass halos are not merging very 
often. We will call this scale m^e/ in the following, our choices for it are 
discussed below. A similar mass dependence was introduced into C("^) by 
FMH05, but without the transition mass rriref- 

With a starburst, /=„ will be larger (more gas will go into stars) for a 
time ~ hurst + i* (that is, after the merger the number of baryons going 
into stars will increase, and these stars will add extra photons for their entire 
lifetime). The average f^Murst at this time is f^^burstt*/ (tburst+t*)- Simulations 
of mergers at low redshifts provide some suggestions for for a starburst, 
e.g. a range of 65% - 85% of the total gas ( Mihos & Hernquist (1994; 1996)) 
or >80% of the cool gas (Springel, Di Matteo and Hernquist (j2005j)) going 
into stars (see also th e studies of gas available at early times by Machacek, 
Bryan & Abel ( 2003[ )). Observationally, starburst star formation is fit by a 



decaying exponential with time scales as short at 20 Myrs in some examples 
but la sting up to hu ndred s of Myrs (e.g.Papovich, Dickinson & Ferguson 
( 2001 ). Shapley et al ( 200l[ ). Conselice |2006l )) with s imilar time ranges used 



to analyze observed samples (Kauffmann et al ()2003| )). 

We only consider Pop II stars, although it would be straightforward to 
generalize beyond this in our formalism. The transition redshift between 
creating Pop HI stars and creating Pop II stars has been estimated to be 

^Note that the photon production rate estimate in FO05 and in FMH05 has C time 
independent, and assumes that all mass going above lO^K is instantaneously c onver ted 



into photons, which is another common approximation, e.g. Somerville & Livio 1)2003 

^"We choose fesc = 0.05 cf. Wyithe & Loeb (,2QQ^)), canonical choices vary by a factor 
of 10. Loeb et al use N^^bary — 4300 for a Scalo IMF of metallicity 1/20 of solar. 
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15 (Yoshida, Bromm & Hernquist ()2004l ). see Fang & Cen ()2004l ) for 



z 

more discussion of constraints on this transition); however a wide range of 
final redshifts where Pop III stars might be important for reionization have 
been suggested (between 2; = 10 and z = 20 e.g. Wyithe & Padmanabhan 
(I2OO6. V Wyithe & Cen (|il06)). (There are even searche s for e vidence of Pop 



III stars at redshifts down to z ~ 6 (Scannapieco et al ((20060).) a result 
we do not consider redshifts before 2; ~ 12 (which are active, for our choices 
of time scales, from mergers at z ~ 15). 

There are expected differences at high redshift, for example the higher 
densities ~ (1 + 2;)^ lead to much shorter cooling times, as well as lower 
metallicities, thus these guesses are just that, however they are a starting 
point for estimating the size of time dependence of 

3.2 Black Hole Accretion 

We now turn to photon production due to black hole accretion of mass Am. 
We take (again assuming a constant phot on pr oduction rate and using the 



notation of Salvaterra, Ferrara & Haardt (|200 



T =^^^''-^)lt. (32) 
which leads to (using the time dependent version of eqn. EH|) 

Ct(mi, mo, ti, m, t) = 6.9 x 10^ ^^^" (33) 

The black hole photon production formula only applies to halos hosting a 
black hole (i.e. with mass above m^/i ,„j„(z)) and includes the effects f rom U V 
photons only. The factor ^ = 0.1 - 0.2 ryd"^ (e.g. Madau et al (|2004l) 1. 



We are assuming that fuy takes into account any local absorption, i.e. it 
is the flux from the whole region due to the mass accretion. The efficiency 
is usually taken to be e = 0.1. The mass accreted Am is often specified in 
terms of halo mass in the range 10~^ — 10~^m. Assuming the mass accretes 
at the Eddington rate 

nibh = rubh/ts , ts = A5Myr{e/0.iy' , (34) 

and commonly used black hole mass to halo mass relations, rribh ~ 10^'^ — 
10~"^m, this gives lifetimes tacc from 0.01 — 10 ts- We will take a relatively 
short accretion time, 0.1 ts, as our fiducial time tacc and rribh ~ 10~^m.^^ 

^^Many different concerns suggest a sfiorter tacc- Relating the black hole mass instead 
to the halo velocity, mhh ^ 10~^(1 + z)'^/^m i.e. including redshift dependence (Bromley, 
Somerville & Fabian l|2004l) .Wvithe & Padmanabhan (,2006i) ). lowers tacc- Accretion at 
super-Eddington rates (argued necessary to get super massive black holes by redshift 6 or 
7, e.g. Haiman l|2004(l ') also gives a shorter tacc- In addition, shorter accretion times are 
also expected from some low re dshift constraints based on the luminosity function (e.g. 
Wyithe & Padmanabhan l|2nn(t V 
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For our purposes we consider all halos above some minimum mass mbh,min{z) 
to host black holes, and all halos below it not to host black holes, and ignore 
the scatter in this relation. We tried to estimate this minimum mass by using 
methods already found in the literature. We have two choices corresponding 
seed black holes at z = 24 in h alos w hich are 3.5a and 3a fluctuations, fol- 
lowing Madau et al's approach (j200J), see the appendix for more discussion 
of this choice and constraints. The two cases used for minimum halo masses 
hosting black holes and m{10'^K) are shown as a function of redshift in Fig. 
1. Our fiducial model takes black hole hosts as 3.5cr fluctuations a.t z = 24, 
concordance cosmology, and the parameters in eq. ESI As mentioned earlier, 
we also explored some models with WMAP cosmological parameters: here 
a detailed application of the black hole halo constraints (some of which are 
described in the appendix) has yet to be done. However, in this case 3a and 
3.5(7 halos at z = 24 have masses 11.3h~^MQ and 600/i~^Mq respectively, 
which is difficult to reconcile (especially in the former case) with initial black 
hole masses of ~ lOOh'^M^. We took one exploratory example: a black hole 
in 5 0" peaks at z = 24 (halos of mass ~ 6 x lO^/i^^M©). The motivation for 
the choice was that at least initially ~ 10~'^ — lO'^^mhaio and this limit 
is also shown in Fig. ^ 



3.3 Recipe for Ionizing Photons 

Combining the above, we have the following recipe for (t{fni,mo,ti,m,t) 
due to a rrii halo merging to mass mo at time which later ends up in 
m at time t < ti + Uon'. Besides the merger being within the Uon of inter- 



Table 1: "Recipe" for (t{mi,mo,ti,m,t) 



cause 


major merger 


mCt 


tion 


quiescent 




re J 


U = 80 Myr 


starburst 


> 1 : 3 


{l3*,sb P*,q) yy^a 

ref 


hurst + U = 100 Myr 


black holes 


> 1 : 10 


Pbhmo 


t,ec = O.U5 = 4.5 Myr 



est (starburst or accretion), the (t{'^i,i^o,ti,m,t) are also zero unless the 
masses also satisfy certain conditions. For the mass range cutoff for a major 
merger we will use 1:3 for starbursts (many ranges are used in th e liter- 
ature) and 1:10 for black hole accretion (e.g. Madau et al (20o3))- For 
black holes mj > 'mhh^min{z{ti)), for quiescent and starburst star formation 
mo > max(m(10*^i^(tj)),mo,mm(^i))- 
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Then (3^^(3bh encapsulate all the non-mass dependent factors: 
3 = f N , ik£j_ 

h'*,q J * ^ ^ -y ,bary j 

_ .36 f* N^,bary fesc 80Myr 
Myr 0.1 4400 0.05 i, 

O _ 1.5 f*,burst N~t,bary fesc thnrsf + ** ( "^^^ 

A^*'*" Myr 0.5 4400 0.05 0.8 ^ ^ 

f^^h =6-9x10^^6^^ 

13.6eV 

_ 15.3 fuv/{hu) / mi,h/m \ / e \ 
Myr O.lrj/d-i V 10"* Aq.I^' 

Our fiducial parameters are shown in the second line of each equation. 
Our fiducial a = 2/3. The reference mass m^e/ is taken to change with 
redshift, 

(e.g. Wyithe & Loeb ( 2003[ )^^). Some numerical studies have found a transi- 



tion m ass rriref which doesn't change significantly with redshift (Keres et al 
diooi), the corresponding transition mass is also almost a factor of 10 smaller 
at z = 0). 

For we estimate assuming an Eddington rate and tacc/ts 1 

so that Am = rribhie^'""'^^^ — 1) ~ rribhtacc/ts and we take the black hole mass 
proportional to its host halo mass, rribh = lQ~^mo. The suggested relation of 
black hole mass to halo velocity mentioned earlier would give a larger ratio 
of black hole mass to halo mass and a bigger effect from black holes. In 
addition, super Eddington accretion will also shorten tacc and increase [3bh- 
Note that changing +tbursti tacc also changes the range of integration for ti 
in calculating not only the prefactors jSbh- For starburst star formation 
f*,burst gives an effective of f*,burstt*/{U + hurst)- 

With these assumptions, the C^tijni^mQ^ti^m^t) in equation El are a set of 
constants for any given mo (or m for quiescent star formation). They are 
only nonzero when tj, m^, mg are in the right range, and for starbursts, black 
hole accretion or their combination. These constants were denoted as Q{mQ) 
in section 12.41 All of these numbers rely on a huge number of estimates of 
unknowns, i.e. the contributions of starbursts, black hole accretion and their 
relative strengths. As mentioned earlier, ideally one could turn this around 
and use it to estimate the contributions from starbursts and quasars, but 
many uncertainties are involved. Here we are interested in the sporadic and 
time dependent nature of the mergers changing photon production rates. 

Using these definitions, Ct,a{jn) (equation ^ is shown in Fig. |21 for two 
black hole assumptions for a series of different times. The black hole and 
star formation contributions are shown separately as well. The quiescent 

-'^^They interpret the m^/'^ behavior as due to a potential weU effect. We have taken their 
default parameter values for = 176 km s"""^ and Acrit = 187r'^+82a;— 39a;^, x = 17m(z) — 1, 

n,niz) = an(0)(l + z)3/[an(0)(l + + Ha]. At Z > 9, Acr^t/18^T^ ^rr^iz) - 1. 
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contributions to (t,a roughly scale with the dotted line in each case, but do 
not change with redshift. The ratio of merger-induced photon contributions 
to total photon contributions for our fiducial model and for the same redshifts 
is shown in Fig. IHl The source mergers have a similar photon production rate 
to the quiescent rate, increasing with increasing mass and redshift. There 
is also a sharp increase in the merger contribution relative to the quiescent 
contribution once the masses are large enough to host black holes. Although 
the black holes have a very large contribution per halo for high mass halos, 
the sharp decline in numbers of halos as mass increases somewhat limits their 
effects. 

4 Consequences for bubbles 

The framework and prescriptions of the previous two sections can now be 
combined to predict the resulting properties of the ionized bubbles. In prac- 
tice there are two steps: C^tijni^mQ^ti^m^t) is used to find the required over- 
densities for given bubble masses M, i.e. 5m- The barriers 5m then determine 
the bubble size distributions and the time evolution of the ionization fraction. 

Our fiducial model, shown unless otherwise stated, is for the 3.5cr seed 
black holes at z=24 as described earlier and in the appendix, with param- 
eters given in eqn. EH We also considered Sa seed black holes ai z = 24, 
a = 0, ruref constant, f^fesc ^f*fesc^ no stars, no black holes. As men- 
tioned above, we also varied the cosmology in several ways. We considered 
ilfe/i^ = 0.0225 leaving all else fixed and cosmological tilt n = 1.05 leaving 
all else fixed. The latter model seemed to have very similar effects to the 
one increasing /*/esc- We took some combinations of the variations with 
each other as well. The two lower integration limits for mg, (^«o.mm = 0.1m 
or mo,mm = Tn{T = W^K)) gave very similar results when a = 2/3. We 
considered minihalo and MHR recombination for all the models, and MIPS 
recombination for the fiducial model and 5 representative variations. One 
other quantity we must fix is an initial condition; this is discussed below. 

4.1 Barriers 

To find the required 5m given a bubble mass M we use the generalization of 
equation 121 equation ITTl 

J dt'iiCf), - A^{1 + 5^Mt'))^C{RUm = 1 • (37) 

Because time dependence is built in, we either need to integrate all the way 
back to the time of the very first ionizing photon production or include some 
initial condition of global ionization fraction. We choose to fix an initial 
condition aX z = 12. At z = 12 halos are only active from mergers since 
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z ~ 15 with our choices of relaxation times. As discussed earlier 2; = 15 is 
the earliest redshift Pop II star dominance (and thus our assumptions) might 
be expected to hold. 

The distribution corresponding to the initial photons is taken to be a time 
independent C('^) similar to the form given by FMH04 & FMH05, so either 
Cfmh is constant or 



^2/3 ui < m. 



Cfmh<x{ '^^^f' " ■■^"'^ (38) 

const m > rriref 

in accordance with the mass dependence of our Q (a = 0,2/3). Then (fmh 
is normalized by setting the global ionization fraction 

lira dmCFMH{m)Pi{m,tbc\M,6M) = Xi^bc ■ (39) 

At a later time, the ions present in a region of mass M with overdensity Sm 
are a combination of the ions produced in that region at the initial time plus 
those produced since, minus recombinations: 

/ dmCFMHim, tbc)Pi(m, tbc\M, Sm) / 
+ IL dt'liCfh - A„(l + 6n,Mit'))^CiRUm = 1 • ^ ' 

For our fiducial model, we start with an ionization fraction Xi{z = 12) of 
almost zero (10~^). This should isolate the evolution of bubble properties 
due to the change in photon production since the initial time. We also 
experimented with some non-negligible initial photon distributions at 2 = 12 
[xi = 10~^, 10~^, 0.5, 0.14) and an initial condition of Xj = 10~^ at z = 16 
when halos are active from mergers since 2; = 23. This last corresponded to 
Xi = 0.16 at 2; = 12. In practice at 2; = 16 Pop II stars are much less likely 
and our assumptions about photon production rates suspect, as a result this 
model was not our fiducial model, but is useful for studying effects of initial 
conditions. 

To find Rion{t) the integral is divided into 700 steps from the initial time 
to the redshift of interest. At each time, the instantaneous change in ionized 
mass fraction {(f)tdt is added to the mass fraction already present (the 
integral up to that time) to find Mion/M, the ionized mass fraction. In the 
case where more photons were produced than hydrogen atoms in the region, 
we took Mion/M = 1 (allowing it greater than one had a neglible effect). We 
took Rion = [i^^(j^f^^Y^\ however replacing ^ Mi^n or even ^ 

had very small effects on the resulting barriers we found; the biggest effect 
is due to the overall coefficient for recombinations, Larger numbers of 

steps gave indistinguishable results. 

We show in Fig. IH examples of the resulting barriers 6m at z=7, 8, 9, 10, 
11 and 11.9 for our fiducial model and the three recombination prescriptions. 
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At large radius, Cmhr{,R) tends to approach a constant, Cmips is constant, 
and Cmh{R) increases exponentially. Thus the barriers for minihalo recom- 
bination are cut off at large radius. MIPS and MHR recombinations allow 
bubbles to get very large, if there is sufficient photon production, to the 
point where consistency questions involving e.g. travel times within the cor- 
responding bubbles arise. For plots of many of the models we thus choose 
to only include the cases with radii < lOOh^^Mpc, models with much larger 
radii do appear but are difficult to interpret. 

Our time dependent approach produces results which interpolate be- 
tween the two limits used in earlier works (e.g. FZH04,FO05,FMH05): they 
included implicit recombinations at early times, and imposed equilibrium 
between instantaneous photon production and recombination rates at late 
times. These two limits, in terms of Sm, give the following: the early time 
approximation produces a barrier similar to ours for small M, approaching 
a constant (horizontal line as function of M) for large M. The late time 
(larger M) limit gives a barrier 6m that is very close to a vertical line (see 
for example FO05, Fig. 7), combining the early and late time limits gives 
the barrier we found which rises at both low and high M. Physically, for 
small bubbles, the photon mean free path exceeds the bubble radius, thus 
photons escape freely and recombinations do not play an important role; this 
is expected to happen at early times. As bubbles grow in size, the barriers 
for ionized regions decrease as a function of bubble radii, since more and 
more mass (and therefore ionizing photons) can be included in the region. 
At late times, bubble radii grow to exceed the photon mean free paths and 
recombinations limit the growth of bubbles. Of the three assumed IGM gas 
density distributions, the minihalo and MIPS models put more stringent con- 
straints on the bubble growth than the MHR model. The stronger effects of 
minihalo recombination compared to MHR recombination were also observed 
in the models equating instantaneous recombination and photon production 
(FMH05 and FO05). 

4.2 Bubble Size Distribution 

The barriers found as a function of bubble radii translate directly into a 
bubble "mass function" nh{M, Sm) which gives the number density of bubbles 
with masses between M and M + dM. Descriptions of how to derive the 
resulting bubble mass function nb{M, 6m) for linear barrier s can be found e.g. 
in Sheth & Tormen and McQuinn et al ((HqoU). As the barriers 

here are non-linear in ct^(M), we find nb{M, 6m) by simulating 4000 random 
walks directly for each (M, 6m) combination. For linear barriers, 4000 steps 
reproduced the ionization fractions to a 5-10% percent once Xj > 0.01. 

The quantity nb{M,6M)dM = Mnf,{M, 6M)d\n M counts the number of 
bubbles with log mass InM. One might want this quantity if one knows 
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one has a collection of bubbles and wants to know the size distributions of 
these bubbles. However, for many questions it is of interest to know the 
fraction of mass in the universe in ionized bubbles with InM in the range 
InM to InM + dlnM, i.e. ^nb{M,6M)dlnM. This peaks at a larger 
mass than the former quantity due to the added factor of M. The peak 
of this distribution corresponds to a mass with a characteristic radius 
Rc = [3Mc/(47rp(H-5ij^Af^(t)))]^/^. This radius (used also in earlier work, e.g. 
FMH05) denotes to the bubble size which contains the largest fraction of the 
ionized mass. The quantity ^nb{M, 6M)dlnM is dimensionless and can also 
be thought of as probability distribution for P{\nM)dlnM = P{lnR)d\nR 
(however it integrates to Xi rather than to one). This distribution and Rc are 
the quantities shown/derived in previous work, e.g. FMH05. The one big 
difference is that they divide by the overall ionization fraction Xi. I.e. our 
plots give the fraction of the total mass which is in the ionized bubbles with 
the given In M. Our motivation for this is that the Monte Carlo calculations 
rely on counting paths crossing the barriers ^M^with more fractional scatter 
in counts when there are fewer counts. By not including their factor 1/xj one 
can read off immediately which curves have the largest such sampling errors 
(i.e. the ones with the lowest heights). 

In Fig. El we show R vs. 3^nft(M, 5jvf) = P(\iaR{M)) for the sequence 
of different redshift barriers of Fig. H] above. The trends shown here were 
seen in our other models as well. For all the models, the radii and ionized 
mass fraction grow with time. The late time shape of the distribution in In R 
has a strong dependence upon the assumed form of recombinations. This is 
also shown in Fig. [7| below, lower left. 

For minihalo recombinations the width of the probability distribution 
for R narrows as R increases and the universe ages. This schematically 
shows the progress of reionization: at early times bubbles are small and have 
a large range in sizes, depending more strongly on the local density and 
ionizing photon production where the bubbles live. As the global ionization 
fraction grows, or as the bubble volume filling factor increases, bubbles grow 
in size, slowed by recombinations. Finally bubbles saturate in radii and have 
similar sizes. In this scenario we expect to find small bubbles with a large 
scatter in size during early reionization, and large bubbles with similar sizes 
at late reionization. The narrowing of bubble width with increased ionization 
fraction was also seen in the studies of limiting cases in earlier work. 

For MHR recombinations, the bubble distributions do not narrow at late 
times, instead recombinations are so weak that the bubble sizes tend to 
have runaway behavior, getting larger and spreading more as reionization 
proceeds. In the examples shown, MIPS recombinations are seen here to 
limit the growth of bubble regions, simply because MIPS recombinations are 
so much larger in number than MHR recombinations. However, in models 
where more photons are present (large Xj) the radial distribution for MIPS 



22 



models moves to large R as is seen in the MHR case above for late times 
(large Xi). 

4.3 Ionization Fraction 

The global ionization fraction is the fraction of mass inside bubbles at any 
given time; time histories are shown in Fig. IHl for the MHR, minihalo and 
MIPS recombination ansatze with the fiducial photon production model. 
Even with negligible photons present initially, the ionization fraction gets 
close to 1 by z = 7 for minihalo and MHR recombinations. For comparison, 
a model similar to FMH05 is also shown (given by equation IHHj) . At early 
times the ionization fraction grows more slowly for this time independent ( 
model as the mergers increase photon production as a function of mass more 
steeply than the m^/^ slope. At late times recombinations become more im- 
portant in the full models and slow reionization, hence, relative to them Xi 
for the FMH05 model increases more quickly. The cases for stars only and 
3.5(7 black holes only are shown as well. At late times the ionization fraction 
in the stars only model increases relative to the black holes only model, in 
part this is due to the rise in the black hole minimum host halo mass and 
the stronger decline in black hole merger photon production with redshift as 
seen in Fig. |21 

4.4 Scatter: model variations 

There is scatter in the results both from model-to-model variations and from 
the merger-induced scatter within one model. The former represents un- 
certainties in the modeling, the latter scatter is expected even if the input 
parameters are perfectly known. 

In order to see similarities and variations between models, we show four 
different comparisons in Fig. [7| The first model to model scatter we consider 
is that due to the unknown initial conditions for our formalism. The two 
top plots show the effects of the ^ = 12 initial conditions on the ionized 
mass fraction for different redshifts. The curves are for the fiducial Xi{z = 
12) = 10~^ model and Xi{z = 12) = 10~^, 10~^. At right we show the fiducial 
model and Xi{z = 12) = 0.5. By z = 10 all the models with Xi at z = 
12 < 0.1 predicted close to identical bubble radii distributions. The model 
with Xi = 0.5 at z = 12 has the radii and ionization fraction changing much 
more slowly between z = 12 and z = 9, and then converging to the other 
cases. In addition, the Xi{z = 16) = 10~^ and Xi(z = 12) = 0.16 models 
(sharing the same Xi at z = 12, not shown above) gave extremely similar 
radial distributions to each other at the redshifts above. 

The bottom two plots consider a sampling of different models: differ- 
ent redshifts, photon production rates, recombinations, cosmologies, etc. At 
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lower left models with similar xt but otherwise widely varying assumptions 
are shown. For low global ionization fraction the different models give dis- 
tributions that are quite similar, irrespective of recombination method (and 
other properties such as redshift or photon production methods). For high 
global ionization fraction the distributions depend much more strongly on 
recombination methods and not only the global ionization fraction Xi. The 
Xi ~ 0.95 models with MHR recombination peak to larger R relative to 
the minihalo recombination models. The characteristic radius Rc for several 
hundreds of our model variants is shown at lower right, the general trend of 
larger radii corresponding to larger ionization fraction is clearly visible, with 
a spread that can be read off of the plot. 

4.5 Scatter: merger induced for a fixed model 

The scatter in the photon production rates for a given model due to mergers 
was described in ^2.41 To propagate the scatter in photon production rates 
calculated from Eq. (j27j) to the scatter in bubble barriers, we replace 

/ dt{Cf)t -> / {Cf)t ± ^A{J{Cf)tdty + al (41) 

in Eq. pn| . We add the scatter of the initial condition {af^, corresponding 
to the first term in equation lin|) in quadrature to A{J{(f)tdtY. We take this 
initial condition scatter to also be Poisson in accord with our assumption (and 
that of FMH05) that the initial condition sources are Poisson distributed in 
the bubbles. The resulting barriers (overdensities) are those required for a 
region of mass M to be ionized if the sources within have the one-sigma 
fluctuation calculated above. 

The corresponding barriers for the fiducial model are shown at z = 10 
in Fig. IHl for all three recombination models. The central line indicates 
the average barriers Sm, while the higher (lower) barrier 6m± corresponds 
to the smaller (larger) photon production rate of one sigma fiuctuations in 
{Cf)t- At small bubble radii, the scatters around the mean values are large, 
while at the largest radii, the three barriers converge to the (recombination 
dependent) limiting bubble radius. There are two effects. Smaller bubbles 
presumably contain fewer halos, so the scatter of photon rates per source 
plays a relatively important role in determining the required overdensity for 
bubbles; conversely, the largest bubbles contain more halos and thus the 
effects of the scatter tend to average out within. In addition the barriers and 
thus radii for the scattered and not scattered cases converge when they hit 
the limiting value set by recombinations. At later times the distribution of 
ionized mass fraction as a function of R becomes more and more skewed to 
large radii, causing the characteristic radius to converge with this limiting 
radius. This trend was seen in all the cases where a limiting radius was found. 



24 



Producing fewer photons can lead to a scatter up or down in bubble radius: 
fewer photons means fewer ions all around, but the bubbles present are often 
required to be larger in order to enclose a sufficient number of sources. As in 
the case without scatter, for low R MHR, MIPS and minihalo recombinations 
are relatively weak and thus give similar results. 

The corresponding radii distributions for these barriers are shown at top 
left and right and bottom left in Fig. IHl At bottom right the scatter is 
shown for the minihalo case for the model with Xi = 10~^ at z = 16, for the 
corresponding fiducial model with Xi = 0.16 at 2; = 12 and for this fiducial 
model with the scatter doubled. The scatter for the model evolved from 
z = 16 converges to that starting at z = 12 hj z = 10. 

The full ionized mass fraction distribution is found by combining the scat- 
tered barriers and the mean one in a distribution, with some weighting. This 
weighting should correspond to outcome of the following procedure. In prin- 
ciple there is a sequence of barriers corresponding to adding and subtracting 
the fiuctuations with a continuous coefficient rather than just the two cases 
above, with larger coefficients having smaller weights in the joint distribution. 
Each random walk used to find the ionized mass fraction then has a differ- 
ent barrier ( "walk barrier" ) sampled from this distribution of barriers. At 
large scales all the walk barriers coincide, but at smaller scales fiuctuations 
between different possible barriers come into play to give a range of walk 
barriers. The first crossings of each walk barrier are dependent both on the 
shapes of the barriers in the distribution and on the scales at which each bar- 
rier appears. A fiuctuation from a higher barrier down to a lower barrier in a 
walk barrier at some mass M will produce first crossings not only of a path 
which was counted in the ionized mass fractions shown above, but also of any 
path which would have had first crossings earlier for this lower barrier but 
didn't for the higher barrier. If there was only one such transition between 
the barriers above and it was at a fixed M for all walk barriers, there would 
be a large spike at the mass scale of the transition, due to all these paths 
being included at the same M. (Likewise a drop would be expected for a sole 
transition from a low barrier to a high barrier in a walk barrier.) We do not 
expect these extreme cases to occur for many walk barriers however, and the 
transition masses M will change as well. We make a rough estimate of the re- 
sulting distribution as follows. Most of the time the walk barrier will sample 
the different barriers often enough that a "pileup" of paths that would have 
crossed earlier, causing such a spike, should be suppressed, similarly for the 
sharp signal associated with a single jump up to a higher barrier. In general, 
we expect each walk barrier to sample our 3 barriers frequently, and effects 
from different walk barriers (different samplings of the barrier distribution) 
to average out, so that the resulting ionized mass fraction distribution be- 
comes close to a weighted average of the corresponding distributions for the 
3 barriers shown above. 
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Some general trends were seen across many models and redshifts. The 
scatter effects for the families of models are summarized in Fig. (as 
mentioned above MHR and MIPS cases with Rc > 100h~^ Mpc are left out 
but models with a = are included). The final scatter (z=7) seems small 
in all other cases and the radii converge to a value independent of initial 
ionization fraction (but dependent upon each model's photon production 
rates). We call the peak (characteristic radius) of the distributions in R 
(such as above) for the scattered model Rscat and the peak for the fiducial 
model Rave in these plots. The scatter decreases with increasing radius and 
decreasing redshift. This is shown in the upper and lower left plots: the 
upper left plot shows the effect of reduced scatter with increasing numbers 
of sources {R) included as mentioned earlier, at lower left is the redshift 
dependence. The ratio Rscat /Rave also decreases as a function of increased 
ionization fraction x,, shown at upper right. As mentioned above, this is in 
part due to the ionization fraction profile becoming limited by recombinations 
as the ionization fraction grows, with the characteristic radius converging to 
the limiting radius. At lower right it is shown that the model-to-model 
scatter is often larger than the merger- induced scatter. The average curves 
for 3^nb(M, Sm), the fraction of total H in ionized H bubbles with a given 
In R are shown for a series of models at z = 9 along with the representative 
scattered curves for one of them. Note that the mass ratios chosen to define 
a major merger will also affect the scatter, for example relaxing the major 
merger mass ratio leads to more mergers overall and thus smaller Poisson 
fluctuations. 



4.6 Bubble Boundaries 

One interesting question is how the ionization fraction decreases from 100% 
as one goes outside a fully ionized bubble at any given time (we thank Evan 
Scannapieco for suggesting this calculation and Steve Furlanetto for discus- 
sions about its interpretation). If the ionization fraction dropped from fully 
ionized to e.g. 90% ionized very slowly in space, or if bubbles were very close 
and connected by high density bridges, relaxing the constraint from fully ion- 
ized to 90% ionized would give a very different distribution of barriers and 
corresponding radii. To calculate the 90% ionized case we relax the complete 
ionization requirement for bubbles to 

/ dm(FMH{m, tbc)Pi{m, tbc\M, 6m) (.^^ 
+ IL dmf)t' - Av{l + 5RMt'))^C{RUm = 0.9 . 

The bubble radii, for minihalo and MHR recombinations, are compared for 
90% ionized and 100% ionized bubbles in Fig. ^2 From Fig. ^2 it is seen 
that these would-be bubbles tend to have very similar sizes compared to 
completely ionized bubbles. The main exception is when the global ionization 
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fraction is close to one for MHR recombination, perhaps indicating that a 
significant amount of bubble merging is occurring and/or that there are large 
regions near the bubbles (perhaps not in 100% ionized bubbles) with high 
average ionization fraction. For the above models and a sample of others, 
the fraction of mass 90% ionized or more is within 6% of the mass in regions 
which are fully ionized Bit z — 7, the difference increases to <15% at z = 11. 

5 Conclusions 

We have calculated contributions to Ct,a{nT') including mergers, a time and 
mass dependent phenomena. We developed a time dependent description 
of the requirements for bubble formation {6m) that generalizes previous for- 
mulations. Using three recombination histories and several different param- 
eterizations of the uncertain photon production, the requirements for bub- 
ble formation, characteristic radii and ionization fraction histories were then 
found. 

Our results can be summarized as follows. The resulting solutions for 
5m from our formalism and calculations interpolate between the two limiting 
cases considered in previous work as expected. In comparison with these 
cases, the increased photon production due to mergers gives a faster rise in 
ionization fraction at early times, and then a slower one as recombinations 
become important. At early times the characteristic radii and distribution 
as a function of bubble radius for ionized mass depends most strongly on the 
global ionization fraction. At later times, the radial distribution for identi- 
cal global ionization fractions depends strongly upon recombination prescrip- 
tions. The width of the distribution for the bubble radii narrows for minihalo 
recombinations at later times just as found in earlier work. For MHR and 
MIPS recombinations the distribution as a function of R widens once the 
ionization fraction is large. In addition, fixing global ionization fraction, the 
characteristic radius at late times tends to be larger for MHR and MIPS 
recombinations relative to minihalo recombinations. We also estimated the 
scatter between models (fairly large for different recombination prescriptions 
as described above, decreasing with time for similar models differing in initial 
conditions) and the scatter of radial distributions of the ionized mass within 
a given model. For the latter, recombinations (except when runaway behav- 
ior is present) limit the size of the largest bubbles independent of whether 
there is scatter present or not. As the global ionization fraction increases, 
the characteristic bubble size tends to this limiting size, resulting in char- 
acteristic bubble sizes with very small scatter as well. The model-to-model 
scatter tends to be larger than the merger induced scatter within a model. 

We also explored variations using the WMAP cosmology parameters, with 
the seed black hole prescription mentioned earlier. Because the estimates of 
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the halo masses containing black holes are even more unconstrained for this 
case, and for simplicity, these models are not included in the plots above. 
Summarizing these cases, using the WMAP cosmology parameters and the 
accompanying family of models reduces strongly the number of photons, 
mostly because high mass halos are much rarer at early times. The WMAP 
models of radii vs. ionization fraction overlapped with those seen in Fig. 
[TJ-i.e. there was no additional spread introduced. The results for recombi- 
nation dependence and ionization fraction followed the trends noted above. 
The WMAP cosmology with 5a seed black holes gave larger scatter simply 
because there were fewer sources and hence the Poisson fluctuations were 
larger. 

We also explored how the ionization fraction changed as the bubbles were 
required to not be fully ionized, but only 90% ionized. These mostly ionized 
bubbles had about the same amount of mass in them as the fully ionized 
bubbles, indicating that slightly relaxing the constraint of full ionization 
doesn't change the mass fraction contained in the relevant bubbles drastically. 

The formalism and tools developed here and results found can be ex- 
tended in future work. Many other parameter choices also fall well within 
the range of reasonable guesses due to the large uncertainties at early times, 
explorations of these would be very interesting. Pop III stars can be included 
in the same framework as well. The observational consequences of the trends 
we have identified in our models so far are also a clear next step. This can 
be done using analytic methods based on the distributions we have found. 
Another route would be to implement the (t nierger prescriptions into nu- 
merical sim ulatio ns of histories or semi-analytic simulations such as those of 
Zahn et al ( 2006[ l. and calculate the observational consequences this way. 

Our methods to include mergers might also be useful to include in other 
descriptions of reionization, e.g. those which include the effects of clustering 
of galaxies such as is done by Babich & Loeb liooa). 

We both thank the Aspen Center for Physics and the January winter 
meeting organizers for the opportunity to present this work when it was near 
completion. We are grateful to many people for helpful suggestions, point- 
ers to literature and explanations, including A. Barth, M. Boylan-Kolchin, 
D.Eisenstein, J. Greene, J. Hennawi, G. Howes, I. Iliev, D. Keres, M. Kuhlen, 
C.P. Ma, A. Meiksin, C. McKee, M. Morales, P. Norberg, E. Quataert, B. 
Robertson, E. Scannapieco, R. Sheth, F. van den Bosch, A. Wetzel and O. 
Zahn. We want to thank S. Furlanetto in particular for numerous explana- 
tions and discussions and both him and E. Scannapieco for comments on the 
draft. We are very grateful to the anonymous referee for many suggestions 
and corrections that clarified the work discussed in this paper. Last but not 
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Appendices 

Extended Press-Schechter definitions 

Two probabilities appear in our formalism in the text. The probability that 
a halo of mass m at time t contains a halo of mass mo at time ti is: 

Pi (mo, ti\m, t)dmQ 
_ [T\ daHmo) \ Sc(ti)-&c{t) r {&c{u)-5,{t)f 1 , (43) 

~ \l 2ti\ dmo I [<j2{mo)-cr2{m)]3/2 ^^t^l 2((t2 (mo)-(72 (m)) J ""'0 

The overdensity for collapse at time ti is taken to be time dependent, Sc{ti) = 
1.68/ D{zi) where D{zi) is the growth factor, and af = a{rniY and cr^ = 
(T(mo)^, the variance of the linear power spectrum smoothed over a region of 
mass m (at z=0). 

The other quantity is the time derivative of the above, the fraction of 
halos that have mass mo that have jumped at time ti from mass mi. 

Pi (mi mo;ti)—dmidti = -^=—^ — ^ 9x^/9 [~ '^^']!^'^ ]|-t-^|— • 
mi ^/27!^ [(yf — o'^yi^ dti dmi mi 

(44) 

The factor of m^/mi is added to convert from the number of points originating 
from mi halos to the number of points in mo halos containing these earlier 
mj halos. 

To write the merger fraction, Eq01 consider the following. The number of 
merged halos is the fraction of recently merged halos of mass mo at ti times 
the number (density) n{mQ,ti) of mo,tj halos. Multiplying by their survival 
probability P2(m, t|mo, tj) of surviving to m,t gives the number (density) 
of m halos at time t which had a merger to mass mo < m at this earlier 
time. Dividing the number of recently merged halos by the total number 
(density), n{m,t) gives the recently merged fraction. We then multiply the 
fraction of halos which has a recent merger by the fraction of the total number 
of m,t halos found in a bubble of mass M,6m, i-e. by nh{m,t\M,6M) = 
-^Pi(m, t\M, 6m) to give a number density. The product of these probabilities 
is then 

Pi{mi ^ mo]ti)—n{mo,ti)P2{m,t\mo,ti)—^-—-^Pi{m,t\M,6M) (45) 
mi n[m,tjm 

which equals Eq. HI because 

Pi{mQ,ti\m,t)n{m,t)—dmQdm = P2{m,t\mQ,ti)n{mQ,ti)^—dmQdm (46) 

P P 

(also this can be used to get an explicit expression for P2). These probabilities 
are being multiplied together assuming that they are independent. There is 
the expectation that overdense regions should have larger merger rates (e.g. 
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Scannapieco & Thacker Furlanetto & Kamionkowski (jioO^ ) , however 



the overdensities here are in fact quite small (the bubble overdensities 5m tend 
to be between and 4 with corresponding physical overdensities 5mD{z)), 
so that this effect is expected to be small. 

Minimum black hole masses 

We describe here how we chose our minimum host halo masses for harboring 
black holes. For rnhh^mini^) niany suggestions of high redshift black hole 
histories are available, see for example the review by Haiman & Quataert 



( 2004 ) and references therein. We take our black holes to be descendants of 
the very first stars, expected to be extremely massive due to the d i fficult y 
of fragmentation (see for example Bromm, Coppi & Larson ((199^ liooi), 
Nakamura & Umemura fjl 999: .20011 Abel, Bryan & Norman ff2000i : i20oi . 
Schneider et al (j2002l )). At the end of the lifetimes of t hese s tars, very 



massive black holes are expected to form, e.g. Madau & Rees (|200lf ). however 
the initial mass distribution of the stars and their resulting black holes are 
unknown (many different exam ples a re considered in e.g . Al varez, Bromm 
& Sh apiro (.2006.1 O'Shea et al ^200^ . Scannapieco et al (jioOQl and Madau 



et al (j2004l) l 



We follow Madau et al (|2004l ) and consider two cases, putting very massive 
black holes in all > So" or > 3.5cr fluctuations at redshift z = 24, correspond- 
ing to masses m^(j,m^,^(j{m > 1.4 x 10^h~^MQ,m > 1.4 x IO^H'^Mq) which 
then grow primarily through mergers. We find the minimum mass for black 
holes at our (later) time of interest by requiring at least 90% of the halos 
with mass mbh.mm(^) to have at least one halo of mass mfyh^miniz = 24) in it, 
using extended Press-Schechter. (In practice it appears that at least 4 or 5 
mergers occur for these high a peaks by 2; ~ 15, we thus require at least the 
corresponding number of paths originally at mass mi to be present in each 
final mass halo of mass m at time t in order for that mass to host a black 
hole. We assume the black holes merge when their halos do, e.g. Mayer e t al 



( 2006| ). there is still discussion on this issue, see e.g. Madau et al (j2004l ) "or 



further discussion including considerations of black hole ejection.) A second 
option for estimating halo masses which can contain seed black holes is using 
the requirement of low angular momentum disks which can collapse, the two 
reference models for minimum black hole mass used by Koushiappas, Bullock 



& Dekel (j200J) at z = 9, 12, 15 are approximately the sa me as t hose we find 



extrapolating the minimum mass halos of Madau et al ( 2004fl for 3a and 
3.5(7 host halos at z = 24. (For some reason Madau et al ( 20041 ) get a much 



lower merger rate th an we do, ours is more in accord with the calculations 
of Koushiappas et al (I2OO4I) .) As we are interested in some reasonable black 



hole mass estimates, the similarity between these two calculations suggests 
these are useful reference masses to consider. This is certainly not the only 
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way to estimate the black hole host masses, other optio ns inc lude tracing the 
luminosity function back in time e.g. Wyithe & Loeb (|2002h . 



Black hole radiation in principle constrains black hole masses and num- 
bers. However, compounding the the black hole/host halo mass relation 
uncertainties and accrection mode/speed uncertainties mentioned earlier are 
even more uncertainties tied to black hole radiation estimates. Harder X- 
ray photons may induce additional ionizat ions ( some estimates have X-rays 



dominate the luminosity, e.g. Madau et al (|200J)) in regions with low ioniza- 
tion fr actions, in addition, efficiencies might be higher (Haiman & Quataert 
( 2004 )). The strongest direct constraints on early black holes is that they 



don't produce photons th at exc eed the unresolved soft X-ray ba ckgrou nd (Di- 
jkstra, Haiman & Loeb (2004), Salvaterra, Ferrara & Haardt ( 2005[ l). This 



translates into a constraint on the ratio $ of power law to multicolor disk lu- 
minosities; although we use Madau et al 's (^OOJ) original model of including 
black holes we take $ < 1 to avoid the constraints their models could not 
satisfy ($ ~ 0.1 is common today). Other black hole halo mass constraints 
requiring more evolutionary assumptions over longer periods of time include 
predicti ng cor rectly black hole to halo mass ratios today (e.g. Ferrarese & 
Merritt ( 2000f)'). the number of intermediate mass black holes today (Yu & 
Tremaine ( 20021 )) and the observed quasar activity observed up to z ~ 6. In 



addition, one expects starbursts and black hole growth can be related using 
the relation between central black hole mass and bulge velocity dispers ion ^\ 
models proposed to do this explicitly include those by Cattaneo et al ( 2005[ l 
and Enoki et al (|2003[ ). 



Other scatter assumptions 

Other assumptions are possible for 

(n(m, t|M, 6M)n{m', t'\M, 6m)) ■ (47) 

and 

{n'^{mo\m,t)n\mo\m,t)) . (48) 

besides those used in the main body of this paper. We list two more possi- 
bilities for each. 

For {n{m,t\M,6M)n{m',t'\M,6M)), one might expect that a large fluc- 
tuation at time t implies a large fluctuation at time t' = t + e, giving a 
correlation between different t, t'. In addition, halos of mass m at time t will 
have a different mass m' at time t', so the fluctuations will be related for 
m ^ m' and t ^ t' . 
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We thank Joe Hennawi for this suggestion. 
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One possibility is to give the scatter a decay time (or characteristic Srei) 
and a decay mass rrirei, 

{n(m, t\M, SM)n{m', t'\M, Sm)) = e-('^(*')-^(*))/'5-'e-('^'("^)-^'("^'»/^'("^-')n(m, t) 

(49) 

Other possibihties are to have the decay depend on t rather than S{t) and m 
rather than a^{m). 

Another possibihty is to say that the relation between the mean values 
of n{m, t\M, Sm) and n(m', t'\M, Sm) 

n{m, t\M, SM)n{m', t'\M, Sm) = 

n(m, t|M, Sm) I dm" Pi{m' , t'\m\ t)^n{m" , t|M, Sm) (50) 
m' < rriji' < t 

suggests a relation between the scatters for m' < m,t' < t as follows: 

{n{m,t\M,SM)n{m',t'\M,SM)) = {n{m,t\M,SM) I dm" Pi{m' ,t'\m" ,t)^n{m" ,t\m,t)) 

= n{m,t\M,SM)n{m',t'\M,SM) 

+y^Q{t - t') Jdm"^Pi{m', t'\m", t)S{m - m")n{m, t\M, Sm) 
= n{m, t\M, SM)n{m', t'\M, Sm) 
+^e{t-t') ^Pi{m',t'\m,t)n{m,t\M,SM) ; 
m' < m,t' < t . 

(51) 

and similarly when {m,t) {m',t') when m' > m,t' > t. This is not 
required: the relation between the average values does not have to be followed 
by the scatter, it is an additional modeling assumption. 

The equal mass and time limit, m m', t ^ t' can be taken in the above. 
Writing a'^{m') — a'^{m) = x, S{t') — S{t) — y we have 

P^i^^y) = -^-k^-'^^^"'^^ - -^\-7^^2^~'^'^"^^]^ ' (52) 
^/2^:x^l^ dm dy^^/2^:x^l^ ^ dm ^ ' 

then m — > m' is a; — > and t — > is y — > 0. The two limits do not commute, 
thus another assumption is involved, we take m ^ m' first, roughly writing 
/ dm = J dm{Q{m — m') + 0(m' — m) + {SoiiTi — m') — 1)) and then recognize 
the Dirac delta function as the a; — > limit of the quantity in square brackets. 
This also assumes taking this limit and taking the derivative commute. We 
then get the full two point function for these assumptions to be: 

(n(m, t\M, SM)n{m' , t'\M, Sm)) = n{m, t\M, SM)n{m', t'\M, Sm) 

- t')^PAm\ t'\m, t)n{m, t\M, Sm) 
+^e{t' - t)^Pi{m, t\m', t')n{m', t'\M, Sm) ^ ' 

-^K'm - m')SD{S{t) - S{t'))n{m', t'\M, Sm) 

The derivative of Soiv) is evaluated by integrating by parts and changing 
variables to t,t'. 
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There is also a range of possibilities for the scatter in Ct,a{^) or n"(mo|m, t). 
Again, we describe two. 

One possibility is to take the same scatter as is used for n(m, t|M, (5m) 
and use it for Pi{mo,ti\m,t) = ^n{mo,ti\m,t) and then define the scatter 

for P\{mi — i> mQ,ti) as the limit. Naive attempts to take this limit give a 
divergence however. 

One also might say that a halo of mass mo in a halo of mass m has either 
had or not had a recent merger, which suggests a binomial distribution, and 
that if the total number of halos of mass mo in m, t halos is that this 
number is then pN, implying a scatter Np — p^N rather than the poisson 
N . Implementing this appears to require yet another assumption because 
the time integral of the rate appears, but it will decrease the scatter from 
Poisson (~ Np). 

None of these scatter prescriptions is required given current knowledge, 
thus in the main body of the text we used the simplest version. It is rea- 
sonable to expect that the trends of largest scatter for smallest radii will 
continue as large numbers of enclosed halos will tend to enclose the mean 
number of mergers. 
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Figure 1: Minimum halo masses as a function of redshift. The lowest solid 
line is the mass of a lO^K halo. The top 3 lines show three models for the 
minimum mass for a halo to contain a black hole. They come from choosing 
black hole hosting halos at z = 24 which are respectively 3a (dotted), 3.5(7 
(dashed), or 5a (dot-dashed, uses WMAP cosmological parameters) peaks; 
see the appendix for details. Our fiducial model uses the 3.5a seeds at z — 24. 
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Figure 2: Redshift and mass dependence of (t,a{fn), lines are z = 
5.9,7.1,8.0,9.0,10.2,11.2,11.9,12.6.13.4, bottom to top (mergers are more 
common in the past, increasing photon production and thus Ct,a{rn)). Top 
left: stars and 3.5 a black holes, top right: stars and 3a black holes. Bottom 
left: 3.5 a black holes only. Bottom right, stars only. The dotted line in each 
case is the slope that would arise from a time independent broken power law 
as described in eq. OH Mergers steepen the slope slightly (m^/^ — >~ m°'^) 
at every redshift. 
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Figure 3: Ratio of photon production rates: merger-produced/ (merger pro- 
duced and quiescent) as a function of mass, for the same redshifts as Fig. 
121 low to high redshift going bottom to top. The small wiggle immediately 
above the sharp rise is a numerical artifact. 
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Figure 4: The barriers 5m found ai z = 7,8,9,10,11,11.9 for our fidicial 
model. The sohd hues (left and right) use minihalo recombination (eqn. 
IT31 in text), the dashes (left) use MHR recombination (eqn. IT^ . and the 
dots (right) use MIPS recombination (eqn. ITB|) . The redshifts for minihalo 
recombination are labelled, for MHR or MIPS they coincide at high z and 
then can be deduced by counting lines down. 
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Figure 5: Top left and right and bottom left: ionized mass fraction ( 
3^nh{M, 5m), i-e. fraction of total mass), in bubbles of radius R for three 
recombination prescriptions. The lines are for different redshifts (11.9, 11, 
10, 9, 8, 7), the smallest height curves are at the earliest time. The MHR 
recombinations are weakest and thus result in the largest bubbles, extending 
at late times to sizes where many other concerns arise as to interpretation 
of the model. At bottom right, the minihalo recombination model is shown 
again (solid line) along with a time independent zeta ~ m^/^ model (eq. 
l38|) having the same Xi at each redshift (dashed line). The time independent 
models have comparable central values of R for high redshifts and low ioniza- 
tion fractions, i.e. at early times when recombinations aren't very important. 
Generally their radial distributions are wider as their recombinations are only 
implicit. 
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Fiffure 6: function of redshift starting with Xi = 10 ^ at z = 12. The 

different hnes of the same type are for MHR (top), minihalo (middle), and 
MIPS (bottom) recombination. The minihalo (center) lines are shown for 
two different Monte Carlo runs to illustrate the run to run scatter. The solid 
line is the fiducial model, the short dashed and dotted lines correspond to 
the models with only black holes and stars respectively. The heavy dotted 
line is the ionization fraction for a C('^) ~ m^^^ model (eq. EHl similar to 
FMH05). 
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Figure 7: Variations due to different modeling assumptions. Top left and 
right: initial condition dependence for minilialo recombination. At left, initial 
(z = 12) Xi = 10-3 (solid), 10-2 (dots), 10-1 (dots), for z = 11.9, 11, 10, 9, 8, 7. 
At right initial z = 12 Xi = 10^^ (solid), 0.5 (dots), for the same redshifts. 
Bottom left: distributions for models with Xi ~ 0.42 (dashed), 0.95 (solid). 
For low Xi the In R distributions are similar, irrespective of recombination or 
other modeling assumptions, for high Xi the profile depends strongly on the 
recombination choice (the rightmost peak is MHR recombination, the peak 
left of it is minihalo recombination, for several models). Bottom right: Rc 
as a function of Xi is shown for all the models to illustrate trends and model 
to model scatter. Open symbols are minihalo recombinations. Squares are 
for n = 1.05 models, f^fesc '^f*fesci VLi^h? = 0.0225, octagons are for the 
fiducial model with varied initial conditions, including an initial condition 
set at z = 16, two variations of mo,mm 3cr or 3.5a seed black holes, a = 2/3 
or a = 0, and stars only and 3.5a black holes only. Filled symbols are MHR 
(octagon/square) recombinations. Filled triangles are MIPS recombinations, 
show for the fiducial model and 4 other cases: the variation taking a = 1, and 
the other 3 combinations of 3.5(T or 3a seed black holes and mo.min = 0.1m 
or m(T = lO^K). 
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Figure 8: Barriers at ^ = 10 with scatter for the fiducial model. The central 
line is the mean, lines above and below are la scatter in the time integral 
of Ct,aijn) as using the replacement of equation UTJ The solid/dashed/dotted 
lines are minihalo, MHR, MIPS recombination respectively. The scatter to 
fewer photons results in a larger bubble size: a larger volume is needed in 
order to enclose a sufficient numbers of sources. 
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Figure 9: Top left to right and bottom left: Radii at z = 10 and their 
one sigma scatter up and down for our fiducial case for MHR (xj = 0.42), 
minihalo (xj = 0.38), and MIPS (xj = 0.08) recombination. Bottom right: 3 
estimates for scatter at z=10 for a model with Xi{z = 12) = 0.16: starting at 
z = 16, starting at z = 12 with usual scatter and starting at z = 12 doubling 
the usual scatter. The lines for the scatter to bigger radii starting at z=16 
and that starting at z=12 are very similar by z = 10. The highest dotted line 
shows how doubling the scatter translates into a change in the characteristic 
radii. Solid lines are the mean, dashed lines are scatter to fewer photons, 
dotted are scatter to more photons. Lines not seen indicate an ionization 
fraction too small to appear. 
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Figure 10: Top left: The ratio Rscat/Rave amongst several models as a func- 
tion of Rave- Top right: The range of the mean value of Rc as a function of Xj. 
Bottom left: Rscat/Rave a fuuctiou of redshift. For all three of these, mod- 
els with Rc > 100 are not shown. Filled symbols are MHR (octagon/square) 
and MIPS (triangle) recombinations, open symbols are minihalo recombina- 
tions, specific models are as in Fig. [7| Bottom right: many different models 
are shown at z = 9. Left to right (smallest to biggest when there is overlap): 
fixed niref, ^th"^ = 0.0225, the fiducial model, n = 1.05,/*/esc — > f*fesc * 4 
(with scatter), 3a black holes, a = 0, 3.5cr black holes, a = 0, 3cr black holes. 
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Figure 11: Mass fraction in bubbles as function of radius for bubbles which 
are 100% ionized (sohd) and 90% ionized (dashed) at redshifts 11.9, 11, 10, 9, 
8, 7 (lowest to highest curves) for minihalo and MHR recombinations (left and 
right). The effects of changing between fully ionized to 90% ionized bubbles 
does not change the distribution as a function of radius that much except at 
large ionization fraction. Note the radial scale for MHR recombinations is 
much larger. 



